package grid;

public class Grid3DBlur {
	public static void blur(float[] values, int nX, int nY, int nZ, int radius) {
		int nYZ = nY * nZ;
		int kernelSize;
		float[] kernel;
		int i, index;
		float sigma = radius / 3f;
		kernelSize = 1 + radius * 2;
		kernel = new float[kernelSize];
		float factor = (float) (1f / sigma * Math.sqrt(2 * Math.PI));
		float sigmaSqr = sigma * sigma * 2f;
		for (i = 0; i < kernelSize; i++) {
			int distance = radius - i;
			kernel[i] = (float) (factor * Math.exp(-distance * distance
					/ sigmaSqr));
		}
		float temp[] = new float[values.length];

		// xdir
		for (int x = 0; x < nX; x++) {
			for (int y = 0; y < nY; y++) {
				for (int z = 0; z < nZ; z++) {
					float sumV = 0;
					float sumKernel = 0;
					for (i = 0; i < kernel.length; i++) {
						int cX = x + i - radius;
						if (cX >= 0 && cX < nX) {
							index = getIndex(cX, y, z, nYZ, nZ);
							sumKernel += kernel[i];
							sumV += values[index] * kernel[i];
						}
					}
					temp[getIndex(x, y, z, nYZ, nZ)] = sumV / sumKernel;
				}
			}
		}
		// ydir
		for (int x = 0; x < nX; x++) {
			for (int y = 0; y < nY; y++) {
				for (int z = 0; z < nZ; z++) {
					float sumV = 0;
					float sumKernel = 0;
					for (i = 0; i < kernel.length; i++) {
						int cY = y + i - radius;
						if (cY >= 0 && cY < nY) {
							index = getIndex(x, cY, z, nYZ, nZ);
							sumKernel += kernel[i];
							sumV += temp[index] * kernel[i];
						}
					}
					values[getIndex(x, y, z, nYZ, nZ)] = sumV / sumKernel;
				}
			}
		}
		// zdir
		for (int x = 0; x < nX; x++) {
			for (int y = 0; y < nY; y++) {
				for (int z = 0; z < nZ; z++) {
					float sumV = 0;
					float sumKernel = 0;
					for (i = 0; i < kernel.length; i++) {
						int cZ = z + i - radius;
						if (cZ >= 0 && cZ < nZ) {
							index = getIndex(x, y, cZ, nYZ, nZ);
							sumKernel += kernel[i];
							sumV += values[index] * kernel[i];
						}
					}
					temp[getIndex(x, y, z, nYZ, nZ)] = sumV / sumKernel;
				}
			}
		}
		values = temp;
	}

	public static int getIndex(int x, int y, int z, int nYZ, int nZ) {
		return x * nYZ + y * nZ + z;
	}

}
